{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "638ca9b8",
   "metadata": {},
   "source": [
    "# 方差分析ANOVA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "637f376a",
   "metadata": {},
   "source": [
    "## single factor ANOVA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb5625b4",
   "metadata": {},
   "source": [
    "<table align='center'>\n",
    "    <tr align='center'>\n",
    "        <th>parameters</th>\n",
    "        <th>explanation</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td align='center'>$SST$</td>\n",
    "        <td align='center'>总体平方和</td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td align='center'>$SSA$</td>\n",
    "        <td align='center'>因素A的平方和</td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td align='center'>$SSE$</td>\n",
    "        <td align='center'>误差平方和</td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td align='center'>$a$</td>\n",
    "        <td align='center'>因素A的水平</td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td align='center'>$n_i$</td>\n",
    "        <td align='center'>因素A第i个水平的样本数量</td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td align='center'>$y_{ij}$</td>\n",
    "        <td align='center'>因素A第i个水平第j个实验观察值</td>\n",
    "    </tr>\n",
    "<table>\n",
    "<style>\n",
    "    td{\n",
    "        align='center'\n",
    "    }\n",
    "</style>\n",
    "计算公式如下：\n",
    "$$\n",
    "    \\begin{array}{ll}\n",
    "    SST = \\displaystyle\\sum_{i=1}^{a}\\sum_{j=1}^{n_i}{{(y_{ij}-\\bar y)}^2} \\\\\n",
    "    SSA = \\displaystyle\\sum_{i=1}^{a}{n_i{(\\bar y_{i}-\\bar y)}^2} \\\\\n",
    "    SSE = \\displaystyle\\sum_{i=1}^{a}\\sum_{j=1}^{n_i}{{(y_{ij}-\\bar y_i)}^2} \\\\\n",
    "    SST = SSA + SSE\n",
    "    \\end{array}\n",
    "$$\n",
    "假设检验统计量\n",
    "$$\n",
    "    H_0 : \\mu_1 = \\mu_2 = \\cdots = \\mu_a \\\\\n",
    "    H_1 : \\mu_i 不完全相等 \\\\\n",
    "    F = \\frac{SSA/(a-1)}{SSE/(n-a)}=\\frac{MSA}{MSE}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cff3f14f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>C</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>164</td>\n",
       "      <td>185.0</td>\n",
       "      <td>187.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>190</td>\n",
       "      <td>197.0</td>\n",
       "      <td>212.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>203</td>\n",
       "      <td>201.0</td>\n",
       "      <td>215.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>205</td>\n",
       "      <td>231.0</td>\n",
       "      <td>220.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>206</td>\n",
       "      <td>0.0</td>\n",
       "      <td>248.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>214</td>\n",
       "      <td>0.0</td>\n",
       "      <td>265.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>228</td>\n",
       "      <td>0.0</td>\n",
       "      <td>281.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>257</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     A      B      C\n",
       "0  164  185.0  187.0\n",
       "1  190  197.0  212.0\n",
       "2  203  201.0  215.0\n",
       "3  205  231.0  220.0\n",
       "4  206    0.0  248.0\n",
       "5  214    0.0  265.0\n",
       "6  228    0.0  281.0\n",
       "7  257    0.0    0.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.stats.anova as ssa\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "data_a = np.array([164, 190, 203, 205, 206, 214, 228, 257])\n",
    "data_b = np.array([185, 197, 201, 231])\n",
    "data_c = np.array([187, 212, 215, 220, 248, 265, 281])\n",
    "geniepig_weight_df = pd.concat([pd.DataFrame({\"A\":data_a}),\n",
    "                                pd.DataFrame({\"B\":data_b}),\n",
    "                                pd.DataFrame({\"C\":data_c})],axis=1)\n",
    "geniepig_weight_df = geniepig_weight_df.fillna(0)\n",
    "geniepig_weight_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73fb474e",
   "metadata": {},
   "source": [
    "## F analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "773bce22",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.002307593325694528"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAICCAYAAAD7xXuXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAxOAAAMTgF/d4wjAABT8ElEQVR4nO3deXxU5d3//9dsyWSZJGRjSxAUUFErbkRFgbpb2loolZZC3W3LXfF7t4233367cGOxtpZfW2hta9VWsdpq0drWApJUISqg1qUKIkFAwiYhgezbzJzfHyczSWCyn5nJHN7PxyOPySRnzvlkjHlzLee6HA0NDQYiIiJR5ox3ASIicmJQ4IiISEwocEREJCYUOCIiEhMKHBERiQkFjoiIxIQ73gUADBs2jLy8vEGdIxAI4HK5LKooehKlTkicWhOlTkicWlWn9RKl1sHWeejQIY4ePRrxe0MicPLy8ti7d++gzlFRUUFhYaFFFUVPotQJiVNrotQJiVOr6rReotQ62DpHjx7d7ffUpSYiIjGhwBERkZhQ4IiISEwMiTEcEZGhKBgMYhjWLDdpGAaBQMCSc0VTb3U6HA6czoG1VRQ4IiLHaG1tZc+ePbS1tVl2Tr/fT0NDg2Xni5a+1OnxeBgzZgxJSUn9OrcCR0TkGHv27MHn85GTk4PD4bDknG1tbXg8HkvOFU291WkYBlVVVezZs4fx48f369wKHBGRToLBIG1tbeTk5OB2W/cn0k734eTk5FBdXU0wGOxX95omDYiIdBIas7GqZWNHofemv+NbChwREYkJBY6IiMSExnBERBLE17/+dd577z1cLhcvvfRSvMvpNwWOiEgPbn30dT6qahz0eQzDOG5c6KScVB664YI+n+Mf//gHFRUVg64lXtSlJiIyxD3xxBPMmDGDyspKZsyYwZIlS+Jd0oCohSMi0oP+tEB60tra2u8bJUPmzZvHvHnzGDt2bEJ2pYWohSMiIjGhwBERkZhQ4IiISEwocEREEsTu3bvjXcKgKHBERCQmbBk4Ta0Bbvz9a7y+uzrepYiISDtbToveeqCGlz6oZFhqEheMzY53OSIigk1bOIdqWwB4bZdaOCIiQ4U9A6fODJx9R5vYe2TwS1KIiMjg2TRwmsOfaxxHRGRosGXgVLa3cABe23UkjpWIiETH2LFju/3efffdx9tvvx1+vn79ei688EIuv/xyPv3pT9PQ0NDjuWtraznvvPMsn4Zty0kDh+paSHY7yU1P5rVdVfEuR0QS2RNfhCO7Bn0at2HAsbuIDhsH8/406HMf6+677+7y/Be/+AXPP/88OTk5LFiwgGeffZb58+dHfG1lZSVf+cpXqK62vnfIli2cQ7Ut5GckUzQumw8rGzhc39L7i0REhii/389tt93G1KlTueCCC3jllVcA+PnPf860adOYPHkylZWV4eNvvPHGLot8PvPMM+Tk5ABw8OBBRo0a1e21PB4Pq1at4qSTTrL85+hTC2fx4sWsW7eO/Px8HnzwQfLy8o475o033uC73/1u+Hl5eTk/+tGPuP76662rto8q61soHJbCpFEZPPPWPnYdbiA3PTnmdYiIDVjUAvEPYrXohx56CJfLxSuvvMKWLVtYtWoVYIbDhg0bWLhwISUlJXzpS1/q8TyrVq3C4XBw2WWX8eCDD/LEE090+f6dd97JzJkzB1xnb3oNnLVr17J582bKysrYsGEDS5YsYcWKFccdd/7557NmzRoAampq+NznPsdnP/vZiOdcvnx5l3PU1dUNelOhUPMvEDSoqm9hUl4yjpZ6AD7YvZ8Rrp77LGMlGs3UaEmUWhOlTkicWk/kOg3DwO/309bWRiAQsOy8fr9/wK99++23ufDCC2ltbWXChAkUFxfzyCOPMG/ePFpbW8nJyaGhoYHW1lYAgsEgbW1t4ecA7777LkuWLOGf//wnra2t3Hjjjdx4443d1mkYBq2trV3OERIMBvH7/ezbt++4TeV60mvglJSUMHfuXJxOJ9OnT6e4uLjXky5btow77rgDr9cb8fuLFi1i0aJF4ecTJ06ksLCwz0V3p7CwkEO1zQQNOGn4MCaeNByoAK/PkvNbZSjV0ptEqTVR6oTEqfVErTMQCNDQ0IDH48Hlcll67oG2HM4++2xee+01brzxRrZv387nP/95AIYNGwaAy+XC7XaHz+90OvF4POHn27dv5+abb+bpp59m9OjRfarT4XCQlJQUseZAIIDb7Wb06NH9eo96HcOpr6+noKAAAIfDQWNjz/e1VFdX89JLLzFr1qw+F2Gl0D04eenJ5KSZ3WhV9ccntIhIorj11ltpaWnhkksu4aabbuLXv/51v15/0003cfToUW6++WYuueQSfve730Wp0p712sLx+XxdQiZS86qzxx57jLlz5/armWWl0JTo/IxkctPNZK5q0KQBEUlcHo+Hhx9+uMvXOk9ZXrx4cZfv/eEPf+jyPDTJoD+isbNory2coqIi1q9fD5gTAdLT03s8/oknnmDOnDnWVDcAoZs+831ehqW1B45aOCIicddr4MycOZONGzdSXFzMggULWLhwIStXrmTlypXHHbtr1y6cTifDhw+PSrF9EWrh5PmS8bicZKV6NC1aRGQI6LVLzev1UlpayurVq5kzZw5FRUXdHjtu3Dg2bdpkaYH9FRrDyfeZ4zc5aUlUNaiFIyISb326DyclJYXZs2dHuxZLHKptweGA7PbutJz0ZMo/rotzVSIiYruVBirrW8hJS8btMn+03PQkjjS24Q8E41yZiMiJzXaBc7i+JTw7DQhPja5uVLeaiEg82S5wmtsCpCZ13IiUk66ZaiJiP9FaLfrAgQNcdtllXHPNNRQVFbFlyxbLarZd4LQFDJLcHT9WTrpu/hSRE8vdd9/N5MmTw89Dq0WXlpYybNgwnn322W5f+8ADD7BkyRLWrFnD7bffzs9+9jPL6rLd9gSt/iAeV0fg5Kbp5k8RGbg7Su+gom5waz2CuTbZsTfEF/oKWXH58WtTHsvv9/P1r3+drVu30trays9//nPAXC36mWeeoba2lnXr1oUXVg6tkzZjxgzAXC06pLfVou+5557wcja9Hdtf9gucQJAkl1o4ImIfsV4tGmDnzp089thjvPrqq5b9HLYKHMMwaAt0beHkaHkbERmEvrRA+qJ1ENsTvPvuu1xyySUAnHHGGZx22mk88sgj4dWe8/PzaWnp+W/cf/7zH+655x7WrVsHwO23387tt98esc7a2lq++MUv8vDDD4f30bGCrcZwAkEDwwCPu3OXmlo4IpLYzjzzzHBLY/v27eHxmbS0tD69fvv27SxYsICnnnoq4n5mnTU0NHDddddx1113hUPOKrYKnLaAAYDH1dFPmpHixu10cFiBIyIJKparRf/4xz/mrbfe4uc//zmXXHIJCxcuHGz5YbbqUmttv7kzuVMLx+FwkJOepC41EUlYsVwtesmSJdx33339LbFPbNXCafWbgdN5DAcgOy1ZC3iKiMSZrQKnLRA5cIalejja2BaPkkREpN0JETiZKR7qW/wEg0Y8yhKRBBK6V8Yw9PeiO6H3pr8bbdpqDCcUOEmurm9ChteDYUBds5/MVE88ShORBOF0OvF4PFRVVZGTk2PZ7sXBYJBAIGDJuaKptzoNw6CqqgqPx4PT2b82i60Cp9Vvpm7npW2AcMjUNrcpcESkV2PGjGHPnj1UV1dbdk6/34/bPfT/5PalTo/Hw5gxY/p97qH/0/dDd11qGV7zx6xpaqMw5lWJSKJJSkpi/PjxBINBy7rW9u3bx+jRoy05VzT1VqfD4eh3yybEVoHT2sMYDpiBIyLSVwP9wxqJw+HA5XL1fmCcRbNOe00aCE2LPqZLLaM9cGoVOCIicWOrwGntbtKAWjgiInFnq8DpWNomcpdabbMCR0QkXmwWOO0tnGO71Lxq4YiIxJstA0eTBkREhh5bBU6LPzSGc+ykAXMyXm2TP+Y1iYiIyVaB010LJ9ntwutxqoUjIhJH9gqc8GrRxy9FkZni0aQBEZE4slfghGapuY//sTK8HrVwRETiyFaB03EfzvE/VmaKRzd+iojEka0Cp7tp0RAKHL+WHBcRiRNbBs6xkwbAXG2gNRCkuS0Y67JERASbBU5rL5MGQKsNiIjEi60CJzRpINIYTuctCkREJPZsFTjdbU8AWjFaRCTebBU4oftwups0AGrhiIjEi70Cpw8tHAWOiEh82CxwQtsT9DBpQIEjIhIXtgqcFn8Qj8uBw3F84HRsUaAFPEVE4sFWgdMWCEbsTgPITNW0aBGReDphAic0LVpdaiIi8WG7wIk0Qw0gLcmNwwF1zepSExGJB1sFTmvAiHjTJ4DT6SA92U1di1o4IiLxYKvAaWufNNCdDK9HLRwRkTixVeC09jCGA+DzuhU4IiJxYqvA6WnSAJiBo0kDIiLxYa/A8Qcj7vYZ4lOXmohI3NgqcFoDBsm9tHDMPXECMaxKRETAZoHTFgjicXc/acDXfi+OWjkiIrFnv8DpsYVjrjZQp9UGRERizlaB0+rvfdIAqIUjIhIPtgmcoGHgD3Z/4yd0buEocEREYs02geMPdr81QUhGuIWjLjURkVizT+C074XT3VpqoC41EZF4sk3gtIVbOL13qWmLAhGR2OtT4CxevJipU6cya9YsKisrezx227ZtXHrppbS0tFhSYF917PapFo6IyFDUa+CsXbuWzZs3U1ZWxp133smSJUu6Pba1tZXbb7+d5cuXk5ycbGmhvWnrU5eaJg2IiMSLu7cDSkpKmDt3Lk6nk+nTp1NcXNztsaEwWr9+PampqZx66qkRj1u+fDkrVqwIP6+rq6OioqK/tXdxuPoIAE0N9d2eq77FXGHgYNXRQV9voKqrq+Ny3YFIlFoTpU5InFpVp/USpdZo1tlr4NTX11NQUACAw+GgsbEx4nE7d+7k8ccf57nnniMjI4P58+fzyCOPRAydRYsWsWjRovDziRMnUlhYONCfwbx+VTPwMTnDMrs9VzBo4HBsIehKHvT1BiOe1+6vRKk1UeqExKlVdVovUWqNVp29dqn5fL4uIdPa2hrxuHfeeYcrr7ySs88+m3HjxjFjxgxeeeUV6yrtRVswCPTcpeZ0OkhP0iZsIiLx0GvgFBUVsX79egDKy8tJT0+PeNxpp53Gjh07CAQC+P1+XnvtNU4++WRrq+1BeFp0D5MGQHviiIjES6+BM3PmTDZu3EhxcTELFixg4cKFrFy5kpUrV3Y57vTTT+dTn/oUl112GUVFRZx99tlMnz49aoUfqy/TokFbFIiIxEuvYzher5fS0lJWr17NnDlzKCoq6vbY4uLiHicVRJO/D9OiwWzh7K5qiEVJIiLSSa+BA5CSksLs2bOjXcugtAZ6X9oG2nf9VAtHRCTmbLPSQGgttZ4mDYDZpdbqD9Li1yZsIiKxZJvAaevHpAHQzZ8iIrFmn8Dpx6QBUOCIiMSabQInPGmgly61jBRtUSAiEg+2CZy2PuyHA2rhiIjEi30CJ9C+0kAvXWrahE1EJD5sFDh9vw8H0NRoEZEYs03g9GdaNKhLTUQk1mwTOP1t4ahLTUQktmwTOOEWjqZFi4gMSbYJnHALx9370jagFo6ISKzZJnBCLRyXs+fASU9y43CohSMiEmu2CRzDzBucjp4DJ7wJmwJHRCSm7BM47Y89x43J3IRNXWoiIrFkn8Bpb+I4emnhgDZhExGJB9sETkhfWzi68VNEJLZsEzh9HcMBdamJiMSDbQInGPqkD00cn9dDiz9Iqz/Y+8EiImIJ2wROxxhO78fqXhwRkdizTeCE9G0MR6sNiIjEmm0Cp79jOKDAERGJJdsETmg0pi9datoTR0Qk9mwTOKE7Px196FQLdalparSISOzYJnAM+j9poFYtHBGRmLFP4IRaOH2cFg0awxERiSX7BU6futQ0hiMiEmv2CZz2x/7dh6MWjohIrNgocNrHcPpwbEeXmlo4IiKxYp/A6cd9OOnJauGIiMSa7QKnL11qLqeD9GRtwiYiEkv2CZz2x77shwNaMVpEJNZsFzh9ZQaOWjgiIrFim8DBMHD2rXEDmBMHtNKAiEjs2CZwgkbfu9NAXWoiIrFmm8Ax6NuU6BBtwiYiElu2CRyMvs1QC9FqAyIisWWbwDEw+rSsTYhWGxARiS3bBE6wny2cDC3gKSISU7YJHAN1qYmIDGW2CRyMvq0UHdKxJ45aOCIisWCbwOl3CydZC3iKiMSSbQInaBh9WrgzRJMGRERiyzaBA/2/DwcUOCIisWKbwDH6eeenJg2IiMSWfQKH/rVwNC1aRCS27BM4hoGzH6t3podaOC1q4YiIxIJ9Aof+tXC0CZuISGzZJ3D6uVo0mOM4ug9HRCQ27BM49K+FA9qiQEQklmwTOAyoheNRl5qISIzYJnCCGP1aaQDUwhERiSXbBI5hDKRLzUNzW5C2gDZhExGJNtsEDvRvLTXQ8jYiIrHk7stBixcvZt26deTn5/Pggw+Sl5cX8bgbbriBvXv34vF4SEtLY9WqVZYW2xPDoF9rqUHX1Qay05KiUZaIiLTrNXDWrl3L5s2bKSsrY8OGDSxZsoQVK1ZEPHb79u28+uqr/R68t0LQMPrdpabVBkREYqfXwCkpKWHu3Lk4nU6mT59OcXFxxON27tzJ7t27ueqqq2hpaeH2229n/vz5EY9dvnx5l9Cqq6ujoqJigD+Cye/3EwjQr/P4m+oA+HDPfjKDtYO6fl9VV1fH5DpWSJRaE6VOSJxaVaf1EqXWaNbZa+DU19dTUFAAmNOOGxsbIx6XnJzMU089xaWXXkp1dTVTp05l1qxZpKWlHXfsokWLWLRoUfj5xIkTKSwsHOjPAIDLXY47SL/OM+awA9iPN2MYhYUjBnX9/hjszxpLiVJrotQJiVOr6rReotQarTp7nTTg8/m6hExra2vE40aMGEFRUREA2dnZZGVlsWfPHovK7J250kD/XtOxCZu61EREoq3XwCkqKmL9+vUAlJeXk56eHvG4hx56iGXLlgGwdetWDh48yNixY62rtBf93YANtEWBiEgs9dqlNnPmTJYtW0ZxcTFlZWUsXLiQlStXArBgwYLwcV/5yle4+eabOf/880lLS+PRRx8lJSUlepVH0P9p0WrhiIjESq+B4/V6KS0tZfXq1cyZMyfcbXaslJQUnnzyScsL7KuB3fipFo6ISKz06T6clJQUZs+eHe1aBsWg/2upaVq0iEjs2GalgYFMGkjXSgMiIjFjn8Ch/zd+upwO0pJc1KpLTUQk6uwTOAPYngDMiQPahE1EJPrsEzj0f9IAaIsCEZFYsU/gDGDxTggFjlo4IiLRZp/AGcAGbBDa9VMtHBGRaLNP4BgDe53P69YmbCIiMWCbwIGBTxoATY0WEYk22wRO0ADnALrUMrTagIhITNgmcMyVBvr/Om0zLSISG7YJHAwDxwAmRoe61HTzp4hIdNkmcNTCEREZ2mwTOMFBrDQAChwRkWizTeDAwFcaAE0aEBGJNtsEzkBWiwZ1qYmIxIp9AmcAq0VD5z1x1MIREYkm+wTOINZSA7VwRESizVaBM5AutfRkBY6ISCzYJ3BgQPfhuF1OUrUJm4hI1NkocIyBTVNDWxSIiMSCfQJn4HmjLQpERGLAVoEzkEkDoBaOiEgs2CdwGNikAQi1cBQ4IiLRpMDBbOE0tQW0CZuISBTZJnAGulo0dOyJU69WjohI1NgmcIIDvA8HOq82oMAREYkW2wSO2aU28EkDoD1xRESiyTaBM4jbcLRFgYhIDNgmcAyMQU0aAC3gKSISTbYJnOCg7sNRC0dEJNpsEzgw8C61DI3hiIhEnW0CZ6CrRQNkpJgtnJomBY6ISLTYJ3CAgbZxMtsDp7ZJXWoiItFin8AxDJwDbOFkqoUjIhJ19gkcBt6llprkwuV0KHBERKLIPoFjDGwDNjBvGM1M8WjSgIhIFNkncBh4CwfMbrVatXBERKLGFoFjGOaUgYHehwPmTDV1qYmIRI9NAqf9k0G0cDK8brVwRESiyB6B0/44iLwhM8VDQ6v2xBERiRZ7BE57E2egq0VD53tx1MoREYkGWwROsL2JM5gWTmi1gVqtpyYiEhW2CByD0KSBgZ9DN3+KiESXPQIn1MKxoEtNgSMiEh22CJyQQXWpeTWGIyISTbYIHCumRauFIyISXbYInKAFN34qcEREossWgWPFfTgZKdqETUQkmuwROOH7cAZ+Dt2HIyISXfYInPbHga4WDeDzqktNRCSa7BE47avROAfx07icDnzJbu36KSISJfYIHEtGcbRitIhINPUpcBYvXszUqVOZNWsWlZWVvR7/u9/9jttvv33QxfVVx42fgzuPNmETEYmeXgNn7dq1bN68mbKyMu68806WLFnS4/E7duzghz/8oWUF9oU17RszcNTCERGJDndvB5SUlDB37lycTifTp0+nuLi422P9fj//9V//xV133cU777zT7XHLly9nxYoV4ed1dXVUVFT0s/QORxrNcZfGhoZBncdDG7VNbXy0Z8+g7unpSXV1dVTOGw2JUmui1AmJU6vqtF6i1BrNOnsNnPr6egoKCgBzrbLGxsZuj/3JT37Cpz71Kc4666weA2fRokUsWrQo/HzixIkUFhb2p+4uvHUtwFZ8vvRBnWdEdjXBnbUMyx8ZXuomGgZTY6wlSq2JUickTq2q03qJUmu06uy1S83n83UJmdbW1ojH/fvf/+b111/vEiSxEpo0YEWXGuheHBGRaOg1cIqKili/fj0A5eXlpKenRzzumWeeoaqqimuvvZa77rqLdevWxW4sx4LVogGyUpMAONqowBERsVqvXWozZ85k2bJlFBcXU1ZWxsKFC1m5ciUACxYsCB+3dOnS8OcbNmzg8ccf57vf/W4USj5e0MJZaqCbP0VEoqHXwPF6vZSWlrJ69WrmzJlDUVFRryedNm0a06ZNs6TAvujoUhtsC8cMHLVwRESs12vgAKSkpDB79uxo1zJgVt2Hk5XS3qXWFHmcSkREBs4mKw2YBjtpQC0cEZHosUXgBNsHcZzOwUWOxnBERKLHFoETYl0LR11qIiJWs0XgWLHFNEB6shuX06EuNRGRKLBH4Fg0S83hcJCV4uGoutRERCxni8Cx6j4cgMxUDzVq4YiIWM4WgRPaYnqQcwYA2ls4GsMREbGaPQKn/XGwXWpgLm+jMRwREevZI3As7FLLSvHQ4g/S3BYY/MlERCTMJoFjzWrRYI7hgG7+FBGxmj0Cp/1xsKtFg5a3ERGJFnsEjpVdamrhiIhEhT0Cx6L7cECBIyISLbYInGDQfLTkPpzwemrqUhMRsZItAifUwrHkPhzt+ikiEhX2CByLtpgGc1o0oOVtREQsZovAsZLGcEREosMWgWPlLDWf14PDoTEcERGr2SJwguG11AafOC6ngwyvRy0cERGL2SJwrNpiOiQr1cMRBY6IiKXsETihpW0sShxzAU91qYmIWMkegdP+aMWNnwDZqR6qGxQ4IiJWskfgWNzCGZaWRIs/SFOrVowWEbGKTQLHfLTiPhyAYe03f1arW01ExDL2CJz2R6smDWSnmYFzRN1qIiKWsUfgWHgfDnTc/HlELRwREcvYInCCFm7ABpAd6lJTC0dExDK2CJxQC8dpxeqdaAFPEZFosEfgYHELJ00tHBERq9kicDrtMW3J6YalaQxHRMRqtgicYGjSgEXny0ppn6WmLjUREcvYInA6NmCzJnKS3E58yW5NixYRsZA9AsfiadEAWWkedamJiFjIHoHT/mhh3pCdmqQWjoiIhWwROEGL11IDc2q0xnBERKxji8DB4rXUwJwa3dQW0AKeIiIWsUXgWH0fDnQs4KlxHBERa9gjcKLQwhmm9dRERCxli8Cx+j4cMPfEATjSoHEcEREr2CJwQhuwOS38adSlJiJiLXsETvujVVtMg5a3ERGxmj0CJwo3fmoBTxERa9kkcIzeD+onBY6IiLXsETjtj1atpQYdYzhV9QocEREr2CNwotCl5nE5yUr1UNXQYt1JRUROYPYInPCNn1ZOjDa71dTCERGxhj0CJwotHIDctGSN4YiIWMQWgRNevNPi82anJVHd2EogaP2kBBGRE40tAifEyqVtAHLSkzAMOKp7cUREBs0WgROtLrWc9qnRVepWExEZNHsEThRWiwbISU8GNDVaRMQKtgicYNB8tLpLLTvcwtHUaBGRwXJbebJAIMCrr74KwMUXX4zL5bLy9N3quPHT2vPmpGu1ARERq/SphbN48WKmTp3KrFmzqKys7Pa4G2+8kb///e8888wzzJs3z7Iie2NEYYtpgJw0s0vtsLrUREQGrdcWztq1a9m8eTNlZWVs2LCBJUuWsGLFiuOOq6yspKioiG984xsAjBs3zvpquxGN1aKhcwtHXWoiIoPVa+CUlJQwd+5cnE4n06dPp7i4OOJxeXl5fOMb3+DgwYP84Q9/YObMmd2ec/ny5V1Cq66ujoqKigGUb6qqqgbg8OFKKiqsC4dA0MABVFTWDKq+zqqrqy05TywkSq2JUickTq2q03qJUms06+w1cOrr6ykoKADMQfnGxsYej3/vvfdYt24d11xzTbfHLFq0iEWLFoWfT5w4kcLCwr7WfJxhBwxgL/n5+RQW5g34PJFkpW6jKegaVH3HsvJc0ZYotSZKnZA4tapO6yVKrdGqs9cxHJ/P1yVkWlt7Hs+44ooreOGFF3jkkUeoqqoafIV90NGlZr2c9GSq6tWlJiIyWL0GTlFREevXrwegvLyc9PT0iMe99957XHvttQSDQZqamjAMg9TUVGur7Ua0bvyE9uVtNEtNRGTQeu1SmzlzJsuWLaO4uJiysjIWLlzIypUrAViwYEH4uDPPPJPp06czbdo0vF4vP/nJT0hJSYle5Z10rKVmfeLkpifx+u42/IEgbpctblsSEYmLXgPH6/VSWlrK6tWrmTNnDkVFRd0ee/fdd3P33XdbWmBfROs+HDBbOIYBRxrbyPMlW38BEZETRJ9u/ExJSWH27NnRrmXgwn1q1p86dC9OVUOLAkdEZBBs0UcUrftwAHLbQ+ZwncZxREQGwxaBEwxGZ6UBgLz00GoDmqkmIjIYtgicjjEc6xMnz2euNlBZp8ARERkMewROFKdF56V7AahUC0dEZFDsETjtj9G48TO3vYVzWC0cEZFBsUfgRGm1aIDUJDdpSS61cEREBskmgWM+Wr0BW0ieL1ljOCIig2SPwInSFtMhChwRkcGzR+BEuYWTm55MdWMr/kAwKucXETkR2CJwgtFbaAAwWziGoa2mRUQGwxaBE+5Si1LihG7+PKRuNRGRAbNH4LS3cKJx4yd0Wt5GM9VERAbMFoETbaEWjiYOiIgMnC0CJ5r34QDhVaJ1L46IyMDZInA6Jg1EuUtNK0aLiAyYLQInPIYTpZ8mN719AU+1cEREBswegUP0tpgGSHa7yEzxUFnXHJXzi4icCOwROFFcLTokNz1JkwZERAbBJoET3aVtAPJ9XgWOiMgg2CNw2h+jtbQNwPCMZGqb/TS1BqJ2DRERO7NH4MSgS214hrkR28e1GscRERkIewROlFeLBgWOiMhg2SJwglFeLRo6BY7GcUREBsQWgdOxllr0rjE8o30BT7VwREQGxB6BE+X7cKCjhXOwRoEjIjIQtggcYjBpIL+9haMuNRGRgbFF4ARDfWpRlOx2MSzVo0kDIiID5I53AVboWEttAE2clnrY+xocfA+GjYXCKeAbEfHQ4RlejeGIiAyQPQKn/bHfcbOjBFbdBk3VHV9zuuGC22D6XZCa3eXw/Awvr++qxjCMqM6IExGxI1t0qQ3oxs/XfgePz4FgAK5cAjf8Az67AkZOhs2/ht9cCh9v6fKS4b5kmtoC1Db7LatdROREYYvACRr9nKV24B1YczfknAJfXQ9T74Rxl8K5X4Fb1sFnfgH1H8PDV8Pul8MvG5FpzlRTt5qISP/ZInBC+tTCaWuCZ24HHDDnEcge1/X7TiecdyMseNZ8/uSX4MB/ALNLDeDjWs1UExHpL1sETr+2mH7td1C5DT75f2Hk2d0fN+5S+NIT4G+GP86Bmn0Mb9/5UzPVRET6zx6B0/7Ya5dawA+bfwvpw+Gib/R+4nHTYNZvze61p29keLoLgI+1EZuISL/ZInCCfW3hbPs71O6FC24Fd3LfTn7mbDOc9r7G+HfuB+BjrTYgItJv9pgWHZql1tuBGx8AVzKcd1P/LnDFYtj7Bmlv/pZLXbnsrxk+gCpFRE5stmjhhLrUnD01cQ6Xmzd4njUH0vP6dwGXB2b9Bjxp3O95kJojhwdcq4jIicoegdOX+3A+WG0+TvrcwC6SPQ6u/iEjqGTu0YcGdg4RkROYTQKnD/fhbF8D7hRz9tlAnXcT21MmMytYQsvu1wZ+HhGRE5BNAqf9k+7yprEa9myCUz4JnpSBX8jh4KXx/0MAJzz/TXOVAhER6RN7BE77KE63a3fuKAUjABOvGfS1kkaczsOBa0mufBfeeGTQ5xMROVHYI3B622J6+xrzceLVg77WyKwUlvtn0+gdAf+6B+orB31OEZETgS0CJ9jbtOiPXoX8M7rddqA/RmWm0IiX9eP+DzTXwLrvD/qcIiInAlsETniL6UiJU7sf6vZDwfmWXGtklrmeWplnKpz8SXjnCdj3piXnFhGxM1sETuhGnIj34YTCYPS5llwqJy2JJLeTAzXNcPW94HCarZwY7DoqIpLIbBE4Pf6p3x8KnPMsuZbD4WBkptcMnOGTYPI82F0G5S9Ycn4REbuyR+D0tJbavn+b99/knW7Z9UZmetl/tMl88sn/Z55/3ffNxUFFRCQiWwROx6SBYxInGIR9b8GoyeCybtm4UVkp1Db7aWjxQ8YouPgb5pYHb//RsmuIiNiNLQKnYy21Y75R/SG01FjWnRYyKtO8efRATXsr5+JFkJoLL94LrQ2WXktExC7sETjhLrVjEmffv81HiyYMhIRmqu072r5NgTcDZtwN9Qdh068tvZaIiF3YI3DaH48bwvn4PfNxRA87ew5AwbBUAPYdaer44nk3wrCx8OpyaDpq6fVEROzAHoHT3aSBqg/B6YZhJ1l6vYJhZpdaxZHGji+6PDDj/5o3g278paXXExGxgz4FzuLFi5k6dSqzZs2isrL7pVy+//3vc8UVVzBt2jR+9atfWVZkb7pd2qZqBwwbZ4aBhUZnmYGzt3MLB+CsL0DuqWa3WoP2zBER6azXwFm7di2bN2+mrKyMO++8kyVLlkQ87q233qKxsZGSkhL++c9/snTpUgKB2KymbBgRutMCfqjeBTnjLb+e1+Mi35dMRXVj1284XfDJ70BrPbz8M8uvKyKSyHqdK1xSUsLcuXNxOp1Mnz6d4uLiiMedc845nHPOOQAcOXKErKwsXC5XxGOXL1/OihUrws/r6uqoqKgYSP0ANDY14nDQ5Rzu2j2MDLZRm5RPzSDO3Z38NBd7DtcfX7fvHIZnn4b7td9x8KRZBNK6bkddXV1teS3Rkii1JkqdkDi1qk7rJUqt0ayz18Cpr6+noKAAMLusGhsbezzeMAy++c1v8v3vd7+o5aJFi1i0aFH4+cSJEyksLOxrzcdJ9n6Mg9qu59j+PgAZY88hYxDn7s7Jww/z3sH95I0YhddzTLBeswSeuJ5RO56AT/9/x712MD9rrCVKrYlSJyROrarTeolSa7Tq7LVLzefzdQmZ1tbWHo+/9957KSws5Itf/OLgq+sjw+D4PrWqHeZjFLrUAAqzQ+M4EQJ4wlVQMAXefBSO7I7K9UVEEk2vgVNUVMT69esBKC8vJz09vdtjH3jgAbZs2cL9999vXYV9YuA8NnGiHDihqdEVx04cAHO63GXfhaAf1v8kKtcXEUk0vQbOzJkz2bhxI8XFxSxYsICFCxeycuVKVq5c2eW48vJy7rrrLioqKrj66qu54oor2L9/f9QK78wwIk2J3gFJ6ZbsgRNJaGr03mMnDoScPB3GTYN3noTK7VGpQUQkkfQ6huP1eiktLWX16tXMmTOHoqKiiMdNmDCB+vp6ywvsi2CkrQGqdkDOKd2s6Dl4he0tnOOmRnd22ffg4SvhpR/BF34flTpERBJFn+7DSUlJYfbs2d2GTbwZHJMrrQ1Quy9q3WlgLm/jcPQSOIVTYMLVsOUZOPhe1GoREUkENllpgK5jOKGB+uxTonbNZLeL4T5v19UGIvnkd8zHl34UtVpERBKBPQIHus5Sq20fO8ocHdXrFmanHH/z57FGTYbTPwPb/qGtqEXkhGaPwDGMrnPU6g6Yj76RUb1u4bBUjjS2Udfc1vOBM74DOODFpVGtR0RkKLNJ4BwzhlMbm8A5KScNgI+qemnlDJ8EZ82BHSUkfaxWjoicmOwROBhdF+6MUQtnbK45U213VR82XZt+NzicZP5bK0mLyInJHoFz7OKddQfA6YHUnKhed2xfWzgAuePh7Hl4D7wGO9dHtS4RkaHIFoETjDSG4xsBzuj+eKHA2X24j9tKT78Lw+mGf/2wY08FEZEThC0CJ+IYTpS70wAyUz1kpXr61sIBGHYSDRNnw97XYEdJdIsTERli7BE4dAqcQBs0VEZtSZtjnZST1rcxnHa1k28HVzL86x61ckTkhGKLwMEAR6hTrf5j8wsZo2Jy6bE5qRyqa6Gx1d+n4wNpI+CCW+DAO+a9OSIiJwhbBE6XMZzwlOjYtXCgjxMHQi75b/Ckwov3QjAYpcpERIYWWwROl5UGwlOiY9fCgX5MHABIz4cpt8OhreY6ayIiJwB7BI5h4DwucGLTwhmb2z5TrT8tHICpd0KSz1xjLdC37jgRkURmj8Ch0xhOKHBiNobTz6nRIanZcNF/mdsovPtUFCoTERla7BE4nadFx3gMZ1iqhwyvm139mKkWdtFC8GbBS/eBv+etu0VEEp1NAqfT9OK6A2ZXVbIvJtd2OByckp/OzsoBbD7nzYSpi+DoR/D249YXJyIyhNgjcOi0tE1olYEYOiUvncP1rRxtHEArZcpXITUXNvwU2pqtL05EZIiwR+AYdEwaaKyCtLyYXn98fjoAHw6klZOcbk6Trt0H/9Y21CJiX/YIHNoHcQwDmo5CSlZMr39KXnvgHBrAOA6YN4L6RpqtnOZaCysTERk6bBE4wWB7l1prPRgBcyA+hk7JM2eqDaiFA+BJgRn/FxoPw6srLKxMRGTosEXghMdwmo6aX4hxC2dMdioel2PggQMw+cuQeyps/CXUHbSuOBGRIcIegWMY5rTo5qPmF2LcwnG7nIzNSWPHoUEEjssNVyyGtkZzmrSIiM3YInCg/T6cOLVwwBzH2VPdSIs/MPCTnHotFF4Ibz4Gh8utK05EZAiwReCYi3c6oOmI+YUYt3AATslPI2j0cxHPYzkccOUScxyq9H+tK05EZAiwReCEVxoIdamlDIt5DaGZaoPqVgMYUwSnfRre/ztUvGZBZSIiQ4M9Aof4ThoAmJBvrmzwwcG6wZ/s8h+AwwXrvq9N2kTENuwROHGeNAAwYXg6TodFgZM3Ec5dAHs2wgerB38+EZEhwCaB075adBxbOF6Pi7E5aXzwsQWBA+Z9OZ5Us5UTaLPmnCIicWSPwAGzTy2OLRyAU0f42F3VQFPrIGaqhfhGmHvmVJXD6w8N/nwiInFmj8AxDPMHaTpqtgrcSXGp47QRGRgGlB+yqJVz8SLIKDA3aWuosuacIiJxYo/AodMstTi1bsBs4QBss2IcByApFa78X2iugReXWnNOEZE4sUXgBEMzueKwcGdnp42wcKZayJmfN28G/ffv4eMt1p1XRCTGbBE45n04jri3cMZkp5LicVkbOA4HXHuf+UOuuVvTpEUkYdkmcJzEZ2uCzpxOBxOHp1vXpRYy6hxzcc9dG2Db89aeW0QkRmwROACpNMdla4JjnTrCx+H6FirrWqw98eXfh6R0eOG74Lf43CIiMWCLwAkaBulG+5IycWzhAJwxKhOALftrrD2xbzhM+zYc2QWbfm3tuUVEYsAWgWMY4KN9t804rKPW2ZmjQ4EThZ07L1wIw8bC+p9AzV7rzy8iEkX2CBwMfKEWTpy71E4f6cPpgHf3WtzCAXAnw7X3Q1sDrP4f688vIhJF9ggcA3xGqIWTFddaUpPcnJKXzntWd6mFTLwKJl0H2/6hddZEJKHYInCCBqQzNFo4AGeNzmTvkSaONrZG5wLX3AdJPvjnXdDaEJ1riIhYzBaBQ+cutTi3cADOaB/HeW9fFMZxADJGwWXfhZo9sP7H0bmGiIjFbBE4hgHpoS61IdLCAaLXrQYw5TYYORk2/korEIhIQrBH4ABpRvvWzt6MuNYCMGmUWcO7+6IYOE4XfObnYAThH/8NwWD0riUiYgFbBE7QMHDjN5+4k+NbDJCe7OaUvDT+s/dodC806hy44Dao2AxvrYzutUREBskWgWMY4KF9kzJXfLYmONY5Y4ZRUd1k/YoDx7rsu5A+wtyore5gdK8lIjIINgkcA7fRvumZ0xPfYtqdMyYLgLf2HInuhbwZMHOZuXDpP76pxT1FZMiyR+DQuYUzRAKn0Fzx4K2Ko9G/2OmfNrcx+OB5eG9V9K8nIjIAtggcDMwxHKenfSe2+Dt1hI/UJFf0Wzgh194Pqbnwz29D/aHYXFNEpB9sEThBw8Bt+IfM+A2Ay+ng7IIs3qmowR+IwQyytByza63pCDyvrjURGXpsETgG7S2cIdKdFnLuSVk0tQX44GOL98fpzhmfg0mfg/f/rq41ERly7BE4BrgJDKkWDnSM47y552jsLjpzGaTlmRMIjlbE7roiIr2wR+DQfh/OEAucc08yA+f1XdWxu2haLlz3K2ipgb9+XTeEisiQYYvACRrgNtqGXJdadloSpw73sWlnFUYsx1QmXg3n3wK7y2DTr2J3XRGRHvQpcBYvXszUqVOZNWsWlZWVPR774x//mKVLl1pSXJ+FZqkNsRYOQNHJ2Ryqa2F3VWNsL3zVDyFnApQugYPvxvbaIiIR9Bo4a9euZfPmzZSVlXHnnXeyZMmSbo/9zne+w9q1ay0tsC8Mht4stZALT84BYPPOqtheOCkVZj9orrX2l1u0jYGIxJ27twNKSkqYO3cuTqeT6dOnU1xc3O2xt99+O6effjp79uzp8ZzLly9nxYoV4ed1dXVUVAx8gDtoGLiMNloDQT4exHmioSDJXOPtX+9VMHWkg+rqGI7nkIfv3DvIeuNnNDy9kOppP+zXq2Nb68AlSp2QOLWqTuslSq3RrLPXwKmvr6egoAAAh8NBY2P3XUNjx46lrKys14suWrSIRYsWhZ9PnDiRwsLCvtQbkWH8Bw8BklLSB3WeaCgExufv4b2Pm8PvY0xrHP19OPoeaeV/Je2Ma2Dyl/r18qH2fnYnUeqExKlVdVovUWqNVp29dqn5fL4uIdPaGqVdLAeh4z6codelBlA0Lpv9Nc1UVDfF/uJOJ8z6DfhGmjeEVn4Q+xpEROhD4BQVFbF+/XoAysvLSU9Pj3pR/WUM0VlqIVPH5wJQtqPnCRdRk5YLn38Y/M3w9I3QGuMJDCIi9CFwZs6cycaNGykuLmbBggUsXLiQlStXsnLl0Nh/JTTd2GUMvRs/Q6aekovTARu2xylwAMZOhU/+Pzi0Ff6+SEvfiEjM9TqG4/V6KS0tZfXq1cyZM4eioqIej1+wYIFlxfVF6O/mUFzaJiQz1cPZhVm8uqMK/6X58Svkkm/Cvjfh3afNzdsu+q/41SIiJ5w+3YeTkpLC7Nmzew2beDAAJ0GcBIdsCwdg2oQ86lr8bD0Ux+6s0HhOzgR44Xuwa0P8ahGRE07CrzRgGAae0PbSQzlwJuYB8NqeGC3k2R1vBnzxCfCkmuM5R3uewi4iYpXEDxwgKRw4Q7NLDeDsgkx8Xnf8Awcgb6LZ0mmsgifnQUt9vCsSkRNAwgdOsHMLZ4hsLx2J2+Vk2sQ8th1q4lBtc7zLMXcJvey78PG7sOoWCAbiXZGI2FzCB45hkBBdagBXTRqOAZS8P0R25Lz023D2l2D7Gnjhu/GuRkRsLuEDB8DjGPpdagAzTs3H5YR1Ww/GuxSTwwGf+QWMuRg2PQCv/S7eFYmIjSV84BhG5zGcod3CyUzxMHlUOq98WEV9iz/e5ZjcyfDFP0L2ybD6Lnj/H/GuSERsKuEDJ2gY5m6fMOQDB+DScRm0+oPxvQn0WKnZ8OW/QGoO/OVm2P1KvCsSERtK+MAx6DyGM7S71AAuGZcBwOr3hki3WkjOKWbouJLgyS/BwffiXZGI2EziB45hJEyXGsBwXxLnjsmiZOvHNLYOkW61kFGTze41fxM8PhtX7dDa6kFEElviBw6JM0st5LNnj6KpLUDpUJmt1tnJ02H276ChkvzVt+jGUBGxTOIHTjBxZqmFzPzEKJwO+Ns7++NdSmRnfA6u+xWu+gPw6GegZl+8KxIRG0j8wMFIqDEcgDxfMhefksv6DyqpaWqLdzmRTZ7HkUt+AEd2m6FTN8TGnEQk4SR+4CTQtOjOPjt5FK2BIP/4zxBt5QANp86Bmcug+kP4w0y1dERkUBI/cABPeFp0YrRwAD511khSk1w89foQH5i/4Fb41E+hagf8/hqo3hnvikQkQSV+4CTIatHHSk9285lPjOKdvTW8f6A23uX0bMptcN0DULMXfv8pbVMtIgOS8IETNDpPGkicwAG4/oJCAP481Fs5AOd82dymuqESfn8t7Pt3vCsSkQST8IFjThpIvC41gHPHZDEhP51n39pHU2sCrNZ85myY+0dobYQ/fBq2r413RSKSQBI+cEig1aKP5XA4mH/hSdQ0tfHsWwkyIH/qNXDD38HtNVck+Pej8a5IRBJEwgdOIt742dnnzyvAl+zm96/swjCMeJfTN4UXwC3rILMA/r4ISpdAMBjvqkRkiEv4wAl2WdomsbrUwJw8MPeCQsoP1fPyjsPxLqfvcsfDrSUw6lwoWwZ/ng8tQ2A3UxEZshI+cBJpA7bu3HDxWJwOeHBDgk05Ts+Hm/4JZ30BPngeHr7KvFFURCSCxA8cOs1SG8JbTPekMDuVmZ8YRVn5Yd6uOBrvcvrHk2KuvXbFYjj0Pjw4A8rXxbsqERmCEj9wjMRb2iaSb3xyPAC//Fd5nCsZAIcDLvlvmPdnc3G7P86B0nsgmAAz70QkZmwQOIm5tM2xTh3h45ozRlDy/iHe21cT73IGZuLV8NUNMHIylP0UHrtOa7CJSJgtAieRdvzsyR2Xm62c+9cm8J38w8bCLS+YS+LsLoMHLtK21SIC2CFwEnC16O6cMSqTWeeMZv32Sl4uT6AZa8dyJ5uLfl6/EjDgz1+Gv90BLfXxrkxE4ijxAyeBl7aJ5FtXTSTJ5eRHq98nGEyQ+3K6M+mz8PVX4eQZ8OZj8OuL4cMX412ViMRJwgdOot+Hc6yCYancNHUsW/bX8uTrNthtM2MUzH8WrvmxuQ7bys/BXxdCY3W8KxORGEv4wAltTxDACU5XvMuxxB2XT2BEhpefrPmAw/Ut8S5n8JxOuPBrsHATnHI5vP1H+NUUeG+V2UQVkRNC4gdO+42fAUfit25C0pPdfP8zk6hpamPp8+/HuxzrDDsJ5q8y79sJBuAvN8MfvwCHE3AquIj0W8IHDu2TBoIOe7RuQq49cwSXnZbPs2/tY+0WG00tdjjgE9fDN16HT3wRdqyDBy6E1Xerm03E5hI+cIIGJDns1cIBcyXpH80+i8wUD9955l17dK11lpYLs39rLgI6cjJs/jWsOBc2PwiBtnhXJ3JCaGxr5M2P32Tl1pXcXXY3z5Y/G9XruaN69hjo6FJL+B/lOMMzvNzzuTNZ9ORbfPOpd/jDjRfgdDriXZa1CqeYofPeX2DdD2B1Mbz2W5jxf+GMWbYZlxOJt5ZACx9Uf8CWqi1sObyFLVVb2Fmzk6DRsdK7z+Pj/NHnR62GhP8rbWDgJmC7Fk7IZ88excvllTz1xl5++eIOFl0+Id4lWc/pNLvZTpsJr66AV38Jq26B9T+BGXfDpM+Zx4hIn7QF2ig/Wh4Ol61VWyk/Uo7f8IePGZU2isvHXM6knEmckXMGk3ImkZmcSUVF9HYgTvzAaV/aJmjDFk7IkuvO5N19tfysZDuTRmZwxaTh8S4pOpLSzICZcjts/CVs/i385SbIvx+m3wWnf1YtHpFj+IN+dtbsDLdatlZt5YPqD2gNtoaPyU/J55KCSzgj5wzzI/cMsr3ZMa814f9KBw2DZPwEnMnxLiVqvB4Xv51/Hp974BUW/ektnvrqRZw5OjPeZUVPajZc/n248L/g1eXw2u/g6Rsh6yS4cCGcMx+S0+NdpUjMBY0gu2t3h1stW6q2sK16G03+pvAxw5KHMWXklC7hkp+aH8eqOyR84HSM4aTFu5SoGpOTyoMLzmPeQ5u56Q+v89RXL2Jcrr1/ZtJy4Mr/hYvvgNcehNcfgjX/Ay/eC+ffCFO+Cpmj412lSFQYhkFFXUWXMZf3q9+noa0hfIwvycfZeWeHg+XMnDMZkTYCh2NojvUmfOCAubSNXcdwOjt/bDa/mDuZ/3riTb78u038+asXUZidGu+yoi8tFz75HXMLhHf+BBt/Ba/8whzrmXgN3pNmwugvqrtNEpZhGOyr32eGS9UWth7eytaqrdS1deyim+ZJC4+3hMZcCn2FQzZcIkn4wAm1cJpsPIbT2bVnjWTZ9WfzzafeYe5vN7Ly1iJOyTtBupc8KXD+TXDuDeb9O689CB/8k7wPnodN95pdbefMh6zCeFcq0i3DMPi48eNwqyX0UdPSsS1JijuF07NPNwMm1wyYkzJOwulI7MkzCf9XOrSWWsMJEjgAs84pAODbT/+H63+zkd/dcD7njhkW56piyOk0996ZeDUc3UPNS78i88O/wfr7YP2PYewlcObnYdJ15niQSJwYhsHBhoNsq97Gpl2bqPjA7CKrbu64yTnZlcyp2ad2jLnknMG4zHG4bNhiT/i/0uG11JwJ/6P0y6xzCkhP9nDHk2/yxQc3cf+cT3Dd5BNwPCNrDLXnfYPMzyyFHSXw9uOw/QVzL55/FsP4y+HMOTDxKvDaeKKFxJ0/6Gd3zW62HdnGtqptbKvexrYj27q0XDxODxOHTeSKMVeEWy4nZ52Mx2n/IQGwQ+C0bzF9IozhHOvKScN56qsXcdtjb3Dnn97mlR2H+cFnziAtOeH/s/afyw2nXmN+NNfAtufh3b9A+TrYvgacHhh3KZz6KTj1WsgsiHfFksCa/E1sP7KdD6o/4P3q99lWtY3yo+W0BDpWBElxpzBx2EROyz6N07JPI7stm0tOu4QkG2yjMlAJ/5fJMAw8joCt78PpyScKsvj7HZfw7af/w1Nv7GXzrmp+Pncy55xIXWzH8mbC5HnmR30lvP83+OCfsGsDfPgv+Oe3YeTZMP5KOOWTUDAF3CfuHwHpXmi8pfxIOeVHy9lWvY0Pqj9gd+3uLnfoZ3uzOW/4eeFwOS37NMb4xnTpFquoqDihwwZsEDihdbeCJ0iTNJJ8n5c/3HgBj27czY9Wb2PObzZyw0VjufPyCWSmnrjvCwDpeXDBLeZHcy18WAofrIbta6Hsp+aHJw3GToWTP2m2gvInacbbCaimpYbyI+XsOLqj4/FoOXWtdV2OK0gv4PIxl3cJl7yUvISaLRYvNggc825aO66l1h9Op4Obpo7jolNy+J+//IdHXtnFs2/t5ZtXTuRLU8bgdiX27BZLeDPM9dnOmGVuj7D/bbPFs/NF87H8BfO45AwouADGXARjLoTR50HSCTD9/ARR01LDrppd7KrZxYdHP6T8aDk7juzgUNOhLsf5knxMyJrAhGETmJA1gfHDxjNx2ER8Sb44VZ74Ev+vdFAtnM5OG5HBswun8te39/HjNdv43nNbePjlXdx66cnMOa8Ar0f/cgfMFkzBeebH9GJoqYePXoU9r8KeTeakgw9L2491m11wo84xV7YeeTbkn26LHWbtKhAMsL9hfzhYdtfuDn/eeYYYmLPETsk6hQtHXcjEYRMZnzWe8VnjyU/NV6vFYgkfOI5Ql9oJOGmgO06ng9nnFnDNmSN4cMNOfv/Kbr771/f42brtfOWisXzh/AJGZaXEu8yhJTndnMk28Srzub/FbAHt2WgGUMVm2PfvjuNdSTD8DDN8Rp4NeadD3qmahh1DQSPI4ebDHDx4kL11e6moq+Cj2o/YVbuLj2o+6rKWGJgrIY/LGseloy9lXOY4xmaOZXzWeArSC2w5BXkoSvjAwW/OCjnRu9QiSU1y83+umMhtl57MU29U8FDZLn5Wsp2fl27n4lNymHNeAVdOGkH6iTirrTfuZBhTZH6AeYdxzV448DYceMf82P827H+r6+vS8s3gyW8PoJwJkD0OMkZrXGgAWgIt7Kvbx956M1BCH3vr9rKvfl+XWWEADhyMSh/FlJFTGJsxlnGZ48IfOd4ctVjiLPH/0rS3cAwFTrfSkt3cNHUc8y88iRe3HWLVm3v517ZDvLKjiiTXuxSdnM0Vpw/nk6fmU5idov8pI3E4zBUMsgrh9M90fL32ABx8Fyq3QeUH5uP+t80uuc5cSZA1htyUkTDqdBg2FoaNM8MoswCST7xxgaARpLq5moMNBznQcICDDQeP+6hsqsTA6PI6t9NNQXoBF4y4gCyyOH3k6RT4Cij0FVLoK8Tr9sbpJ5LeJP5f6fZms8ZweudxObnqjBFcdcYIqhtaef4/+1n3/iE2fVhFWflhfsAWRmZ6mTIumynjshnubiZvREDjPj3JGGl+hLriwGwN1e6DQ9ugagcc2QVHdkP1LrwHNsPesuPPk5wBvpGQMcpsDWWM6vjwjYC0PEjNTYjp24ZhUN9Wz+Gmw8d9HGo8FA6Tjxs/pi14/O6uDhzkpuQyIm0Ek/Mnh4Mk9JGfmh/uAquoqKCwUEsZJYqED5zQGM6JttLAYGWnJbHgorEsuGgs9S1+Xi6v5OUdh3ltVzXPvb2f597eD4Br1YeMz0tn0qgMJo3M4OS8NE7KSaMwO4Vkt4IoIofDbLVkFsCEK7p8a++ejyjMdJsBFAqi2v1mQNXuh4rXoNNqwMfxZnaET1qu+XnoIzUbUrLAO8w8LvQxyJAyDIMmfxNHW45ytOUoNS011LTUcLTlKEeaj3QESvNhqpqqONx0+Liurs58ST5Gpo3kolEXMTJtJCPSRjA8dXiXzz2akGFLCf9X2hHUpIHBSk92c82ZI7nmzJEAVDe08u+PjrDx/Qr2NsDWA7U8+9Y+nn1rX/g1DgeMzPAyJieVUZkp5GUkk+/zku9LJt+XzPAML3m+ZFKTXOqi68zhNLdUyBxt3vtzLMOAltquIVR3EBoOQ0Ol+dFYBdUfwt7XoNPNh90xPKk0ebOoT/FRn5xOQ3I69UleGtzJ1LtcNDid1DugwWFQawSpCbZS2VxHw+YANf4GjrbV0Rb093gNl8NFtjeb3JRcxmWOIzclN/yRk5JDrtf8PC81jzSPzbfVkG71KXAWL17MunXryM/P58EHHyQvLy/icX/7299YunQpGRkZLF26lClTplhabES68dNy2WlJXDlpOKf5WsPdFTWNbbx/sJaPqhrYXdXInqpGdlc1sGVfLZt2Vnd7Lo/LQYbXQ0ZK+4fX3f7oISPFTXqSm5QkF8keFyntH16PkxRPx9e8HifJHhcepwOPy4nb1f7odOByJk6YBY0gbcE2GtoaaAu00Rbs+PAH/ebnnb+e4qUtqZCW7Hya/c20+FtoDjTT5G+iJdBCc1sjTS01NLfU0NJaT1NrvXlcoJnmQCsNwVbqg200EMCMpSbzI1BpftoNh2GQEQySFQwyMhDk9GCQrECAzGCQrKBBFm4yHS6yHB6ynEnkOpPJcntx+ZOguQXqqsHdAO6PzckXbm/HoyvJnGbu8pjLDbnc5mP4axG+1+V556+5cTVUQl2SGeROl/kY/tx1zNcT53fFrnoNnLVr17J582bKysrYsGEDS5YsYcWKFccdd+jQIe666y5efvllnE4ns2bNYv369VEpurN9zXv5R84w9jnepWzT0vDXjx1o7I5hRD4u0uu7O2d354j0+vr6etIrum4nEOn1Ea/fj1q7raGv18KgsaGR1N1db3g0DANSIT0VziyEMwF/IEizP0hzW4DmtgAt/iDNbX6a/UH8AYM2f5CGYJCj/iBtzQb+hmP+Vd7t34G+/VxOh4HD4cDpAIfDgQNwOhw4HObfGAehvzUOwGj/mnnR0N8gBwY4DCAIjmD7tQ1wBDEIhj83vx76Wsf3DIIY7Y+0f24+BgjixyCAQaBPP8/AOHCRhMuRjAsPTkcabnJJdqWQ5kjF7UjB40jBTQpuvCQ7kkgNOkgxXKQakGYYpAUN0oMGvqAfZ1Mt6R6D5GATScFGkoxmkoJNJAWbcBstuIOtuI1WPMFWgkYjtcGjuI3Wjq8brb2XbIFR/Tg2iAMDJ4bDSbD90cBJsP0x0vNjjzN/Uxztj4DDfG60/z4ZOI75mgPDYT4GAgF2u9ydvk6X8+Gg07mPPcZp/kY6On+/6/XCNXX+XnuN5vOQY74feu4ATp5B7tmf7se72j+9Bk5JSQlz587F6XQyffp0iouLIx63efNmLr74YnJzcwHIzs7udkBv+fLlXUKrrq6OioqKAf0AFUcq+FOGD9gJH+wc0DkkykK/ZUnmr3Y026LtMUE40oy+/qu2/X96wwk4MYxOzw1Hx9cMZ/jrBp2e4w6/FsP8AxF6HYYbw3BBpw+jy+fuCF93d3luGEkQ9GAYHvMxmASGByPoAcNjHt99csecgyBJ+EnCTzJtJNNKsqONJPy48eMhgJsAbkcAD37cBDq+hh+PI/S5+X0Xwfbv+9tf0/G99ljA1R4pLoK4HB2fO9s/XGZsdDk+/LkjeMx5guZ5MNo/N2swf7YucdLxucMIf7/9twMnwU5REOE1x8dV9+fvdK7Q15yOvv9jsy/KAkk4Cy+29Jyd9Ro49fX1FBSYK+s6HA4aGxsjHldXV9clXDIzMzlw4EDEwFm0aBGLFi0KP584ceKAZ5rMy/kmn6z8PLV1jYw9ZXyX74X+FdvlaxGa1ZGO68/ruxPp2P379jN69PHbCAy21r5evz/X2rd3X/i/fa81RPzS4N/Xvhw7FGYqhVqOoQakcezX259XVFRQUFAYblVGarR2nKNv5wyfo5vXdT62r/bt2xfx93Sw+tob0Ff79u9n9Kj+tHM61WJpJb3bv28/o0aPCv+jKFyHFYUY7Wc97pehf88nJ3k4euRI1P5/6jVwfD5fl5BpbY3cVM7IyKChoWN2TSAQsPyXK2J9qZn4TvoEFRUV5KbkRv16g1XjrkmYtZg8Ts8Jv7ptX4VC8Pjc7PoFj8tJknvor2vXkOImO23o/7dvTHGTk54c7zL6pCnVTW4C1Hr0yJGonbvX3/yioqLwWEx5eTnp6ZG3Mz7vvPN45ZVXMAyDQCDA66+/zpgxY6ytVkREElavLZyZM2eybNkyiouLKSsrY+HChaxcuRKABQsWhI8bOXIk5513HvPnz8fv9zNlyhRGjhwZvcpFRCSh9Bo4Xq+X0tJSVq9ezZw5cygqKur22F/+8peUlpbS1NTEtddea2mhIiKS2Pp0H05KSgqzZ8/u0wkvv/zyQRUkIiL2NPRHL0VExBYUOCIiEhMKHBERiQkFjoiIxIQCR0REYkKBIyIiMaHAERGRmFDgiIhITChwREQkJhQ4IiISEwocERGJCQWOiIjEhAJHRERiQoEjIiIxocAREZGY6NN+ONF26NAhRo8ePahz1NfXd7v99VCSKHVC4tSaKHVC4tSqOq2XKLUOts7Kyspuv+doaGgwBnzmIWTChAmUl5fHu4xeJUqdkDi1JkqdkDi1qk7rJUqt0axTXWoiIhITChwREYkJ2wTOHXfcEe8S+iRR6oTEqTVR6oTEqVV1Wi9Rao1mnbYZwxERkaHNNi0cEREZ2hQ4IiISEwocERGJiSFx46ed1NfXc9NNN9Hc3MyRI0f40Y9+xKWXXhrx2BtuuIG9e/fi8XhIS0tj1apVMa42MTzyyCM89dRT4edvvPEG7777LiNHjjzuWL2nvWtra+P666/nv//7v5k2bRrPPPMMv/jFL0hOTuakk07iN7/5DS6XK+Jrt2/fzqc+9SnGjx8PwFe/+lVmzZoVy/KHpM7v6YUXXshnP/vZ8PcqKyu5+OKLWbFiRcTXnkjvacIEzuLFi1m3bh35+fk8+OCD5OXlRTzub3/7G0uXLiUjI4OlS5cyZcqUmNb55JNP8oUvfIHrr7+eF198kfvuu6/bwNm+fTuvvvoqDocjpjWGnHXWWeEVHi666CJ+8IMfRDwu3u/pzTffzM033wzA5s2b+e1vfxsxbCB+7+mxf8TLy8u55ZZbcLvdzJ8/P1x/JM3NzcybN4+amhomTZrE8uXLo1Z/a2sr119/Pfv27Qtf+09/+hNr1qwhJSWFK6+8ks2bN3PxxRdHfP2rr77Kt7/9bb72ta9Fpb7Ojn1PH374YX7961+Tm5sLwK9+9StOOeWUiK+N53ualJTEmjVrwt+fM2cOCxcu7Pb1sXpPI/1j+LTTTuPLX/4ygUCAGTNm8L3vfa/b1xuGwde+9jW2bdtGYWEhDz30EF6vt181JESX2tq1a9m8eTNlZWXceeedLFmyJOJxhw4d4q677uL555/nz3/+M8XFxTGuFG677Tauv/56AD7++ONu/zDu3LmT3bt3c9VVVzFt2jQef/zxWJbJhx9+yCc+8QnWrFnDmjVrug2bofCedvaDH/yAe+65J+L34vWetra28oUvfIG9e/eGv3brrbeyZMkSSktLWbVqFRUVFd2+/oc//CHnnHMOpaWl5Ofn8/TTT0e13gceeIBzzjkHAK/Xy1NPPUVKSgrBYJCqqipGjBjR7WvLysp45JFHuOKKK/j85z9PdXV1VGqM9J6+8sorPProo+Hf2e7CBuL7nna2bt06CgsLOf3007t9baze09A/hv/+979zzz33cN999/Gtb32LL3/5y5SWlvLRRx+xefPmbl//0EMPEQgEWL9+PVdeeWW3LbaeJETglJSUMHfuXJxOJ9OnT2fTpk0Rjwv9yyw3N5fs7Gyys7N7/B89mqqrq7n//vu5++67I34/OTmZp556inXr1vHXv/6VpUuX0tDQELP6Xn75Zd544w2uuuoqLr/8ct56662Ixw2l97SkpISJEyd2u+5ePN/Tzn9wjh49yv79+5kxYwYOh4PLL7+cDRs2dPvakpIS5s2bB8C1117Lv/71r6jVmZSUxKhRoyJ+b8WKFVx00UWcfPLJ3b5+9uzZlJaWUlJSwvnnn89vfvObaJV63B/xV199lTvvvJPp06fz7W9/u8fXDpX3dNmyZXzzm9/s8fWxek8j/WP4pZdeCn+tt/fJivc0IQKnvr6egoICABwOB42NjRGPq6uro7CwMPw8MzOTAwcOxKTGztra2rjxxhv5zne+E+6XPdaIESMoKioCIDs7m6ysLPbs2ROzGs8880z+8Y9/8MILL/C///u/fOc734l43FB5TwF+/etfc+utt3b7/Xi9p8f+wamvr+/ynmVlZfX4nnU+Pl7vb2lpKc899xw//elPezzusssuw+fzAWaX7Pvvvx+Veo59T9va2rj33nspKSnhpZde4v3336esrKzb1w+F9/TNN98kKyury+9CJLF6T0M6/2M4OTmZlJQUoPf3qfPfgqysLA4ePNjvaydE4Ph8vi4h09raGvG4jIyMLv+iDQQCGEZs72sNBALcfPPNXHnllXz+85/v9riHHnqIZcuWAbB161YOHjzI2LFjY1QlTJo0iQkTJgDmL/nWrVsjHjcU3lMw/0VWUVHBJz7xiW6Pifd7GuLz+bq8Z36/v8f3rPPxvR0bDZs2beJ73/sef/rTn8J/fCIJBoNceumlVFdXYxgGzzzzDGeffXZMavR4PFx77bWA+Y/OM844o8c/zPF+TwEef/zxcOuhO7F+T4/9x7Db3TGM39v71PlvwUDf04QInKKiItavXw9AeXl5t0tnn3feebzyyisYhkEgEOD1119nzJgxsSyVxx57jOeff57nnnuOK664gi996UusXLmSlStXdjnuK1/5Cv/5z384//zz+frXv86jjz7a4//sVrvtttt4+eWXAVi1ahWTJ0+OeNxQeE8B1qxZw7Rp08LPh+J7GpKZmQkQ/hfgxo0bewy+zr/fmzZtinlIzp8/n4aGBubNm8cVV1zB6tWrAbjmmmu6HOd0OlmyZAlXXnklU6ZMwefzxWy5lq1bt/KFL3yBYDBIXV0dJSUl3f7OQvzfU4DVq1czY8aMLl+L53sa6R/Dp512Wrg7fdOmTZx00kndvn7KlCnh93Tjxo09HtudhFjaprm5mcsuu4ypU6dSVlbGLbfcEp4dsWDBgi7HfuMb3+DIkSP4/X6Sk5N57LHH4lHykPfhhx9y880309jYSEFBAb/4xS/Cv0x6Twfm9ttvZ/78+UybNo0nn3ySBx54gKKiIkpKSnj55ZdJT0/nmmuu6TKDCWDbtm3MnTuX6667jieffJKnn366xz+mJ5LO7+m9997Ln//8Z7xeLzfffDNf/epXAfSe9tHvf/97vvWtb3HuuecCkJeXx9e+9jW+9a1vcfXVV/P000/z4osvMnr0aG644QZ+/OMfd5lAcujQIa6++mpmzpzJc889x/33339cgPYmIQIHoKmpidWrVzN69OhwP313SktLaWpq4tprr+32fgLpH72n/ffuu++yZcsWrrnmGrKysno89sCBA7z00ksUFRX1OGgvfaf3tG927tzJ5s2bmTFjRrezakOOHDnCCy+8wKRJkzjrrLP6fa2ECRwREUlsCTGGIyIiiU+BIyIiMaHAERGRmFDgiIhITChwREQkJv5/O45+yvUI6z4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 480x640 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "\n",
    "# 卡方分布的双侧p检验值\n",
    "fig = plt.figure(figsize=(6,8),dpi=80,\n",
    "                 facecolor='whitesmoke',\n",
    "                 edgecolor='grey')\n",
    "ax = plt.subplot(1,1,1)\n",
    "f_x = np.arange(0,20,0.1)\n",
    "f_y = stats.f.pdf(f_x,3,28)\n",
    "chi_y = stats.chi2.pdf(f_x,3)\n",
    "chi_y2 = stats.chi2.pdf(f_x,28)\n",
    "ax.plot(f_x,f_y,label='f')\n",
    "ax.plot(f_x,chi_y,label='chi2-1')\n",
    "ax.plot(f_x,chi_y2,label='chi2-2')\n",
    "ax.legend()\n",
    "ax.grid(alpha=0.4)\n",
    "p = stats.f.pdf(6.03,3,28)\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5b0f86e",
   "metadata": {},
   "source": [
    "# 单因素各水平均值比较"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be4c04b9",
   "metadata": {},
   "source": [
    "## 因素均值原理\n",
    "通过各因素样本均值和样本方差估计总体均值的置信区间<br>\n",
    "假设各因素以及全体数据误差方差均为$\\sigma^2$<br>\n",
    "那么各因素均服从$N(\\mu_i,\\sigma^2)$的正态分布<br>\n",
    "因素样本均值$\\mu_i$是总体均值$\\mu_i$的一个无偏估计<br>\n",
    "因素的样本方差${s_i}^2$在假设中即为总体方差，可以与SSE产生联系，即$SSE = (n-a){s_i}^2$<br>\n",
    "那么可以使用样本均值构造$Z$，使用SSE构造$\\chi^2$<br>\n",
    "然后使用公式$\\frac{Z}{\\sqrt{\\chi^2/\\nu}}$构造$t(\\nu)$学生检验\n",
    "\n",
    "检验公式\n",
    "$$\n",
    "    \\frac{\\sqrt{n_i}(\\bar y_i-\\mu_i)}{\\sqrt{SSE/(n-a)}} \\sim t(n-a)\n",
    "$$\n",
    "\n",
    "## 因素间差值原理\n",
    "检验公式\n",
    "$$\n",
    "    \\frac{\\bar y_i - \\bar y_j -(\\mu_i -\\mu_j)}{\\sqrt{(\\frac{1}{n_i}+\\frac{1}{n_j})SSE/(n-a)}} \\sim t(n-a)\n",
    "$$\n",
    "Bonferrroni定理<br>\n",
    "设$E_i$为m个随机事件，且$P(E_i)=1-\\alpha$\n",
    "$$\n",
    "\\begin{equation}\n",
    "    P(\\displaystyle\\cap_{i=1}^{m} E_i) = 1 - P(\\cup_{i=1}^{m} \\bar E_i) \\ge 1 - \\sum_{i=1}^{m}P(\\bar E_i) = 1-m\\alpha\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8332dfef",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "288px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
